#include "IonoCorr.h"
#include<complex>

IonoCorr::IonoCorr(EphData data, BLH recPos, SowTime obsTime, Azel satAzel) : data(data), recPos(recPos), obsTime(obsTime), satAzel(satAzel) {}

double IonoCorr::ionmodel() {
    double Pi = consts.get_Pi();
    double CLIGHT = consts.get_SpeedOfLight();
    double tt, f, psi, phi, lam, amp, per, x;
    int week;

    if (recPos.h < -1E3 || satAzel.ele <= 0) return 0.0;

    /* earth centered angle (semi-circle) */
    psi = 0.0137 / (satAzel.ele / Pi + 0.11) - 0.022;

    /* subionospheric latitude/longitude (semi-circle) */
    phi = recPos.b / Pi + psi * cos(satAzel.azi);
    if (phi > 0.416) phi = 0.416;
    else if (phi < -0.416) phi = -0.416;
    lam = recPos.l / Pi + psi * sin(satAzel.azi) / cos(phi * Pi);

    /* geomagnetic latitude (semi-circle) */
    phi += 0.064 * cos((lam - 1.617) * Pi);

    /* local time (s) */
    tt = 43200.0 * lam + time.Sow2gtime(obsTime).time + time.Sow2gtime(obsTime).sec;
    tt -= floor(tt / 86400.0) * 86400.0; /* 0<=tt<86400 */

    /* slant factor */
    f = 1.0 + 16.0 * pow(0.53 - satAzel.ele / Pi, 3.0);

    /* ionospheric delay */
    amp = data.ion[0] + phi * (data.ion[1] + phi * (data.ion[2] + phi * data.ion[3]));
    per = data.ion[4] + phi * (data.ion[5] + phi * (data.ion[6] + phi * data.ion[7]));
    amp = amp < 0.0 ? 0.0 : amp;
    per = per < 72000.0 ? 72000.0 : per;
    x = 2.0 * Pi * (tt - 50400.0) / per;

    return CLIGHT * f * (fabs(x) < 1.57 ? 5E-9 + amp * (1.0 + x * x * (-0.5 + x * x / 24.0)) : 5E-9);
}
